Last updated: 2024-07-17

Checks: 6 1

Knit directory: CocaineIVSA_RNASequencing_mDANeurons/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20240716) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Using absolute paths to the files within your workflowr project makes it difficult for you and others to run your code on a different machine. Change the absolute path(s) below to the suggested relative path(s) to make your code more reproducible.

absolute relative
/home/avm27/Documents/Raw_Sequencing_Data/CIVSA2/CocaineIVSA_RNASequencing_mDANeurons .

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 7a9931d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  data/CIVSA_normalized_counts_DESeq2.txt
    Untracked:  data/DAvsnDA_normalized_counts_DESeq2.txt
    Untracked:  data/gene_count_matrix_CIVSA.csv
    Untracked:  data/sample_info_CIVSA.csv
    Untracked:  data/sample_info_CIVSA_DAvsnDA.csv
    Untracked:  output/GO_all_Combined_Comparisons.txt
    Untracked:  output/GO_all_Combined_Comparisons_Reduced.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_BACKGROUND.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_CRVvsFT.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_CRVvsMN.txt
    Untracked:  output/Gene_Matrix_EntrezIDs_MNvsFT.txt
    Untracked:  output/Sig_Res_Craving_vs_FoodTrained.csv
    Untracked:  output/Sig_Res_Craving_vs_Maintenance.csv
    Untracked:  output/Sig_Res_DA_vs_nDA.csv
    Untracked:  output/Sig_Res_Maintenance_vs_FoodTrained.csv
    Untracked:  plots/
    Untracked:  results/

Unstaged changes:
    Modified:   analysis/_site.yml

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/index.Rmd) and HTML (docs/index.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 7a9931d avm27 2024-07-17 Testing
html bfc97c1 avm27 2024-07-17 Build site.
html 1ee1d08 avm27 2024-07-17 Build site.
Rmd a6b1eb6 avm27 2024-07-17 Testing
html 41cc00d avm27 2024-07-17 Build site.
Rmd 871d786 avm27 2024-07-17 Testing
html 741e5c2 avm27 2024-07-17 Build site.
Rmd 18fb11e avm27 2024-07-17 Testing
html 23fe19a avm27 2024-07-17 Build site.
Rmd 52486bd avm27 2024-07-17 Testing
html 6284345 avm27 2024-07-17 Build site.
Rmd 98eff28 avm27 2024-07-17 Testing
html 5bc128a avm27 2024-07-17 Build site.
Rmd 292e3b5 avm27 2024-07-17 Testing
html 12fbe39 avm27 2024-07-17 Build site.
Rmd 479c9dd avm27 2024-07-17 Testing
html 25a76da avm27 2024-07-17 Build site.
Rmd 89ef098 avm27 2024-07-17 Testing
html cb49269 avm27 2024-07-17 Build site.
Rmd ca44efa avm27 2024-07-17 Testing
html 21b70e3 avm27 2024-07-17 Build site.
Rmd ddb0db7 avm27 2024-07-17 Testing
html a1c090c avm27 2024-07-17 Build site.
Rmd adbd778 avm27 2024-07-17 Testing
html dc04b15 avm27 2024-07-17 Build site.
Rmd 65a97a6 avm27 2024-07-17 Testing
html e58ce5e avm27 2024-07-16 Build site.
Rmd ad8ae96 avm27 2024-07-16 Start workflowr project.

Cocaine craving drives repressive epigenetic differential gene expression in ventral tegmental area dopamine neurons

Tate A. Pollock1,2, Alexander V. Margetts1-3, Samara J. Vilca1-2 & Luis M. Tuesta1-3*

1 Department of Psychiatry & Behavioral Sciences 2 Center for Therapeutic Innovation 3 Sylvester Comprehensive Cancer Center University of Miami Miller School of Medicine, Miami, FL 33136 * Corresponding author ()

Abstract

Dopamine (DA) signaling plays an essential role in reward valence attribution and in encoding the reinforcing properties of natural and artificial rewards. The adaptive responses from midbrain dopamine neurons to artificial rewards such as drugs of abuse are therefore important for understanding the development of substance use disorders. Drug-induced changes in gene expression are one such adaptation that can determine the activity of dopamine signaling in projection regions of the brain reward system. One of the major challenges to obtaining this understanding is the inherent cellular heterogeneity in the brain, where each neuron population can be defined by a distinct transcriptional profile. To bridge this gap, we have adapted a virus-based method for labeling and capture of dopamine nuclei, coupled with nuclear RNA-sequencing and bioinformatics analyses, to study the transcriptional adaptations, specifically, of dopamine neurons in the ventral tegmental area (VTA) during cocaine taking and cocaine craving using a mouse model of cocaine intravenous self-administration (IVSA). Our results show significant changes in gene expression across non-drug operant training, cocaine taking, and cocaine craving, highlighted by an enrichment in expression of repressive epigenetic modifying enzymes during cocaine craving. Immunohistochemical validation further revealed an increase of H3K9me3 deposition in DA neurons during cocaine craving. These results demonstrate that cocaine-induced transcriptional adaptations in dopamine neurons vary by phase of self-administration and suggest that such an approach may be useful in identifying relevant phase-specific molecular targets to alter the behavioral course of substance use disorders.

Differential Gene Expression Analysis with DESeq2

0.1 Setup

Building DESeq2 input files from gene count matrices output by StringTie following removal of samples that are outliers.

Running DeSeq2 analysis with factor levels of FoodTrained vs Maintenance, Craving vs Maintenance and Craving vs Foodtrained

1 Plot Estimates and Gather Results

1.1 Gather Results

These results include overview for each result file and dispersion estimates based on count values.

1.1.1 Craving vs FoodTraining

Craving vs FoodTrained CIVSA Results Overview
baseMean log2FoldChange lfcSE stat pvalue padj
Mrpl15|Mrpl15 2162.26767 0.3962276 0.3468265 1.1424376 0.2532722 0.9784152
Gm39587|Gm39587 20.30129 -2.1136742 3.3911278 -0.6232954 0.5330904 1.0000000
Lypla1|Lypla1 994.90242 0.0198799 0.8799100 0.0225931 0.9819748 1.0000000
Xkr4|Xkr4 4459.42848 -0.0337912 0.2715766 -0.1244260 0.9009780 1.0000000
Gm39585|Gm39585 5797.38480 0.0861245 0.2414588 0.3566843 0.7213282 1.0000000
Tcea1|Tcea1 2657.79689 -0.4037632 0.3330614 -1.2122786 0.2254058 0.9485680
Rgs20|Rgs20 422.79650 -0.1889366 0.7619761 -0.2479561 0.8041684 1.0000000
Gm16041|Gm16041 21.30106 7.5945314 3.6799019 2.0637864 0.0390380 0.6410945
Atp6v1h|Atp6v1h 10659.10495 -0.4956230 0.2967682 -1.6700675 0.0949060 0.8425404
LOC108167788|LOC108167788 11.33585 -0.8220772 3.4983015 -0.2349933 0.8142140 1.0000000
Npbwr1|Npbwr1 924.87540 -0.3079659 0.6636912 -0.4640199 0.6426335 1.0000000
Oprk1|Oprk1 15661.74234 -0.1442457 0.3231485 -0.4463757 0.6553259 1.0000000
St18|St18 279.69289 -2.0082890 0.9250037 -2.1711146 0.0299225 0.5845690
Rb1cc1|Rb1cc1 12217.86115 -0.0796796 0.2844821 -0.2800866 0.7794110 1.0000000
4732440D04Rik|4732440D04Rik 5369.16278 -0.4918180 0.2427412 -2.0261004 0.0427545 0.6632640
Gm26901|Gm26901 163.46252 -0.2293404 0.9562622 -0.2398300 0.8104620 1.0000000
Pcmtd1|Pcmtd1 4321.05905 0.2034726 0.2645104 0.7692425 0.4417494 1.0000000
Rrs1|Rrs1 625.86245 -0.4850563 0.6605362 -0.7343371 0.4627433 1.0000000
Mybl1|Mybl1 184.22714 1.4451732 1.2066675 1.1976566 0.2310507 0.9561403
Adhfe1|Adhfe1 452.68051 -0.1536903 0.7879557 -0.1950494 0.8453543 1.0000000

out of 23922 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 250, 1% LFC < 0 (down) : 48, 0.2% outliers [1] : 2272, 9.5% low counts [2] : 1352, 5.7% (mean count < 2) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.1.2 Maintenance vs FoodTraining

Craving vs FoodTrained CIVSA Results Overview
baseMean log2FoldChange lfcSE stat pvalue padj
Mrpl15|Mrpl15 2162.26767 0.6589546 0.3707556 1.7773286 0.0755142 0.3783101
Gm39587|Gm39587 20.30129 -1.7619611 3.6258386 -0.4859458 0.6270056 0.9347192
Lypla1|Lypla1 994.90242 0.0914200 0.9407122 0.0971817 0.9225821 1.0000000
Xkr4|Xkr4 4459.42848 0.6341912 0.2902594 2.1849118 0.0288953 0.2180911
Gm39585|Gm39585 5797.38480 1.1798031 0.2580070 4.5727554 0.0000048 0.0001973
Tcea1|Tcea1 2657.79689 0.1319203 0.3560044 0.3705581 0.7109667 0.9585345
Rgs20|Rgs20 422.79650 0.3482376 0.8144256 0.4275868 0.6689520 0.9470469
Gm16041|Gm16041 21.30106 7.1396504 3.9105846 1.8257246 0.0678918 0.3581243
Atp6v1h|Atp6v1h 10659.10495 0.1185235 0.3172387 0.3736100 0.7086945 0.9577587
LOC108167788|LOC108167788 11.33585 -0.5355535 3.7403736 -0.1431818 0.8861466 0.9975241
Npbwr1|Npbwr1 924.87540 -0.6449281 0.7097813 -0.9086292 0.3635459 0.7812766
Oprk1|Oprk1 15661.74234 0.2927281 0.3454541 0.8473720 0.3967878 0.8053987
St18|St18 279.69289 -1.0473877 0.9882176 -1.0598756 0.2892012 0.7121028
Rb1cc1|Rb1cc1 12217.86115 0.3066224 0.3041177 1.0082359 0.3133412 0.7349733
4732440D04Rik|4732440D04Rik 5369.16278 -0.1463996 0.2594995 -0.5641616 0.5726441 0.9134234
Gm26901|Gm26901 163.46252 -0.7535570 1.0236447 -0.7361509 0.4616389 0.8512966
Pcmtd1|Pcmtd1 4321.05905 0.5461501 0.2827538 1.9315391 0.0534164 0.3104809
Rrs1|Rrs1 625.86245 -0.1126443 0.7061219 -0.1595253 0.8732551 0.9954319
Mybl1|Mybl1 184.22714 -2.0915101 1.2989231 -1.6101878 0.1073569 0.4494496
Adhfe1|Adhfe1 452.68051 -1.0496649 0.8432838 -1.2447350 0.2132292 0.6257432

out of 23922 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 628, 2.6% LFC < 0 (down) : 1046, 4.4% outliers [1] : 2272, 9.5% low counts [2] : 2713, 11% (mean count < 6) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.1.3 Craving vs Maintenance

Craving vs Maintenance CIVSA Results Overview
baseMean log2FoldChange lfcSE stat pvalue padj
Mrpl15|Mrpl15 2162.26767 -0.2627270 0.3466449 -0.7579138 0.4485026 0.8663335
Gm39587|Gm39587 20.30129 -0.3517131 3.3977053 -0.1035149 0.9175543 0.9936441
Lypla1|Lypla1 994.90242 -0.0715401 0.8799229 -0.0813026 0.9352013 0.9960828
Xkr4|Xkr4 4459.42848 -0.6679824 0.2714805 -2.4605170 0.0138737 0.1490430
Gm39585|Gm39585 5797.38480 -1.0936785 0.2412810 -4.5328003 0.0000058 0.0002852
Tcea1|Tcea1 2657.79689 -0.5356835 0.3330620 -1.6083595 0.1077565 0.4660069
Rgs20|Rgs20 422.79650 -0.5371742 0.7618248 -0.7051152 0.4807386 0.8842358
Gm16041|Gm16041 21.30106 0.4548810 3.4906497 0.1303141 0.8963179 0.9904732
Atp6v1h|Atp6v1h 10659.10495 -0.6141465 0.2967694 -2.0694406 0.0385048 0.2712464
LOC108167788|LOC108167788 11.33585 -0.2865238 3.5021032 -0.0818148 0.9347940 0.9959694
Npbwr1|Npbwr1 924.87540 0.3369622 0.6640182 0.5074592 0.6118327 0.9373635
Oprk1|Oprk1 15661.74234 -0.4369738 0.3231413 -1.3522683 0.1762895 0.5981155
St18|St18 279.69289 -0.9609013 0.9259217 -1.0377781 0.2993734 0.7593673
Rb1cc1|Rb1cc1 12217.86115 -0.3863020 0.2844705 -1.3579687 0.1744736 0.5959401
4732440D04Rik|4732440D04Rik 5369.16278 -0.3454183 0.2427896 -1.4227062 0.1548213 0.5640322
Gm26901|Gm26901 163.46252 0.5242166 0.9577924 0.5473176 0.5841606 0.9276800
Pcmtd1|Pcmtd1 4321.05905 -0.3426775 0.2644209 -1.2959545 0.1949912 0.6289246
Rrs1|Rrs1 625.86245 -0.3724120 0.6606706 -0.5636879 0.5729666 0.9245810
Mybl1|Mybl1 184.22714 3.5366833 1.2149357 2.9110045 0.0036027 0.0605854
Adhfe1|Adhfe1 452.68051 0.8959746 0.7889417 1.1356665 0.2560962 0.7114764

out of 23922 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 1073, 4.5% LFC < 0 (down) : 319, 1.3% outliers [1] : 2272, 9.5% low counts [2] : 2933, 12% (mean count < 7) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.2 Annotate Results Files and Pull Significant DEGs

This subsection is for annotation the results files to include all necessary information for downstream analyses.

After annotation I subset the results to only view the significantly differentially expressed genes with an adjusted p-value > 0.05.

## Subset for significant genes only
results_sig <- subset(res, padj < 0.05)
results_sig2 <- subset(res2, padj < 0.05)
results_sig3 <- subset(res3, padj < 0.05)

summary(results_sig)

out of 247 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 215, 87% LFC < 0 (down) : 32, 13% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 2) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

summary(results_sig2)

out of 1301 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 454, 35% LFC < 0 (down) : 847, 65% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 6) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

summary(results_sig3)

out of 1050 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 854, 81% LFC < 0 (down) : 196, 19% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 7) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

1.2.1 Craving vs FoodTrained

Significant DEGs Craving vs FoodTrained
baseMean log2FoldChange lfcSE stat pvalue padj symbol description entrez ensembl
Gm28836|Gm28836 24.57062 18.512895 3.6825760 5.027159 0.0000005 0.0000748 Gm28836 predicted gene 28836 102635593 NA
Vwa3b|Vwa3b 68.09876 9.322079 2.2726944 4.101774 0.0000410 0.0048096 Vwa3b von Willebrand factor A domain containing 3B 70853 ENSMUSG0….
4930439A04Rik|4930439A04Rik 239.71050 6.086391 1.0543294 5.772760 0.0000000 0.0000025 4930439A04Rik RIKEN cDNA 4930439A04 gene 78119 NA
LOC108167622|LOC108167622 55.54376 21.271643 2.7752955 7.664641 0.0000000 0.0000000 LOC108167622 NA NA
Pard3bos2|Pard3bos2 48.12983 -22.863164 3.6323379 -6.294338 0.0000000 0.0000002 Pard3bos2 par-3 family cell polarity regulator beta, opposite strand 2 105243940 ENSMUSG0….
Col4a3|Col4a3 29.96706 19.601203 3.5629545 5.501390 0.0000000 0.0000079 Col4a3 collagen, type IV, alpha 3 12828 ENSMUSG0….
Dnajb3|Dnajb3 70.36681 9.093523 2.0877638 4.355628 0.0000133 0.0017489 Dnajb3 DnaJ heat shock protein family (Hsp40) member B3 15504 ENSMUSG0….
Asb18|Asb18 12.85015 18.822115 3.6835623 5.109759 0.0000003 0.0000497 Asb18 ankyrin repeat and SOCS box-containing 18 208372 ENSMUSG0….
Gm41920|Gm41920 91.39199 20.766926 3.6195797 5.737386 0.0000000 0.0000030 Gm41920 predicted gene, 41920 105246666 NA
LOC105246710|LOC105246710 41.29964 21.372039 3.6789651 5.809253 0.0000000 0.0000022 LOC105246710 NA NA
LOC108167732|LOC108167732 59.78860 21.982685 3.2366754 6.791749 0.0000000 0.0000000 LOC108167732 NA NA
Fam129a|Fam129a 58.26340 9.542099 2.6179113 3.644928 0.0002675 0.0256134 Fam129a NA NA
Fmo1|Fmo1 51.26689 21.615385 2.7780474 7.780783 0.0000000 0.0000000 Fmo1 flavin containing monooxygenase 1 14261 ENSMUSG0….
Gorab|Gorab 666.84329 2.774594 0.5327144 5.208407 0.0000002 0.0000307 Gorab golgin, RAB6-interacting 98376 ENSMUSG0….
Dnah14|Dnah14 41.60299 21.334681 3.4180513 6.241767 0.0000000 0.0000003 Dnah14 dynein, axonemal, heavy chain 14 240960 ENSMUSG0….
Atf3|Atf3 433.56152 -2.808330 0.7372279 -3.809311 0.0001394 0.0145805 Atf3 activating transcription factor 3 11910 ENSMUSG0….
B230208H11Rik|B230208H11Rik 48.28040 21.761939 2.7106546 8.028297 0.0000000 0.0000000 B230208H11Rik RIKEN cDNA B230208H11 gene 320273 ENSMUSG0….
LOC108167797|LOC108167797 82.55280 9.456510 2.0988898 4.505482 0.0000066 0.0009335 LOC108167797 NA NA
Fam26f|Fam26f 45.72461 21.035286 2.7780345 7.572003 0.0000000 0.0000000 Fam26f NA NA
Mical1|Mical1 109.45019 10.585487 2.6621850 3.976240 0.0000700 0.0077657 Mical1 microtubule associated monooxygenase, calponin and LIM domain containing 1 171580 ENSMUSG0….

1.2.2 Maintenance vs FoodTrained

Significant DEGs Maintenance vs FoodTrained
baseMean log2FoldChange lfcSE stat pvalue padj symbol description entrez ensembl
Gm39585|Gm39585 5797.38480 1.179803 0.2580070 4.572755 0.0000048 0.0001973 Gm39585 predicted gene, 39585 105243853 NA
Snhg6|Snhg6 522.81858 1.260214 0.4234933 2.975758 0.0029227 0.0446701 Snhg6 small nucleolar RNA host gene 6 73824 ENSMUSG0….
Trpa1|Trpa1 47.61966 -22.462571 3.6001396 -6.239361 0.0000000 0.0000001 Trpa1 transient receptor potential cation channel, subfamily A, member 1 277328 ENSMUSG0….
Rpl7|Rpl7 2486.96588 1.005240 0.3038678 3.308150 0.0009391 0.0184408 Rpl7 ribosomal protein L7 19989 ENSMUSG0….
Gm28836|Gm28836 24.57062 20.693225 3.9082681 5.294730 0.0000001 0.0000064 Gm28836 predicted gene 28836 102635593 NA
Gm32790|Gm32790 44.07127 9.563380 2.2173604 4.312957 0.0000161 0.0005913 Gm32790 predicted gene, 32790 102635458 NA
Fer1l5|Fer1l5 230.33092 -2.687277 0.8236933 -3.262473 0.0011044 0.0209149 Fer1l5 fer-1-like 5 (C. elegans) 100534273 ENSMUSG0….
Sema4c|Sema4c 1298.41081 -1.822125 0.6102459 -2.985887 0.0028276 0.0434624 Sema4c sema domain, immunoglobulin domain (Ig), transmembrane domain (TM) and short cytoplasmic domain, (semaphorin) 4C 20353 ENSMUSG0….
Vwa3b|Vwa3b 68.09876 8.717777 2.3896076 3.648204 0.0002641 0.0065801 Vwa3b von Willebrand factor A domain containing 3B 70853 ENSMUSG0….
4930439A04Rik|4930439A04Rik 239.71050 7.561383 1.1113671 6.803677 0.0000000 0.0000000 4930439A04Rik RIKEN cDNA 4930439A04 gene 78119 NA
Il1r2|Il1r2 73.49166 9.935852 2.8549805 3.480182 0.0005011 0.0112029 Il1r2 interleukin 1 receptor, type II 16178 ENSMUSG0….
Il1rl2|Il1rl2 153.74171 -24.241136 3.7848523 -6.404777 0.0000000 0.0000000 Il1rl2 interleukin 1 receptor-like 2 107527 ENSMUSG0….
LOC108167622|LOC108167622 55.54376 21.760261 2.9333442 7.418244 0.0000000 0.0000000 LOC108167622 NA NA
Gm10561|Gm10561 31.31746 -8.823060 2.8132289 -3.136275 0.0017111 0.0295377 Gm10561 predicted gene 10561 628004 NA
Clk1|Clk1 11902.12420 1.056339 0.3336618 3.165897 0.0015461 0.0272603 Clk1 CDC-like kinase 1 12747 ENSMUSG0….
Stradb|Stradb 1233.07664 1.338207 0.3402418 3.933107 0.0000839 0.0024506 Stradb STE20-related kinase adaptor beta 227154 ENSMUSG0….
Fam117b|Fam117b 1077.36822 1.740967 0.5934892 2.933444 0.0033522 0.0491723 Fam117b family with sequence similarity 117, member B 72750 ENSMUSG0….
D230017M19Rik|D230017M19Rik 353.30341 1.817985 0.5121487 3.549721 0.0003856 0.0090159 D230017M19Rik RIKEN cDNA D230017M19 gene 320933 NA
Arpc2|Arpc2 3532.07869 -1.115469 0.2391510 -4.664286 0.0000031 0.0001315 Arpc2 actin related protein 2/3 complex, subunit 2 76709 ENSMUSG0….
Cnppd1|Cnppd1 1597.07545 -1.205071 0.3232300 -3.728214 0.0001928 0.0049685 Cnppd1 cyclin Pas1/PHO80 domain containing 1 69171 ENSMUSG0….

1.2.3 Craving vs Maintenance

Significant Craving vs Maintenance
baseMean log2FoldChange lfcSE stat pvalue padj symbol description entrez ensembl
Gm39585|Gm39585 5797.38480 -1.0936785 0.2412810 -4.532800 0.0000058 0.0002852 Gm39585 predicted gene, 39585 105243853 NA
Trpa1|Trpa1 47.61966 22.0664403 3.3933288 6.502889 0.0000000 0.0000000 Trpa1 transient receptor potential cation channel, subfamily A, member 1 277328 ENSMUSG0….
Rpl7|Rpl7 2486.96588 -1.0763939 0.2841797 -3.787723 0.0001520 0.0050633 Rpl7 ribosomal protein L7 19989 ENSMUSG0….
1700001G17Rik|1700001G17Rik 704.98525 1.5933616 0.4255143 3.744555 0.0001807 0.0057918 1700001G17Rik RIKEN cDNA 1700001G17 gene 67503 ENSMUSG0….
Fer1l5|Fer1l5 230.33092 2.7365066 0.7714631 3.547165 0.0003894 0.0107340 Fer1l5 fer-1-like 5 (C. elegans) 100534273 ENSMUSG0….
Sema4c|Sema4c 1298.41081 1.9826581 0.5709146 3.472776 0.0005151 0.0133359 Sema4c sema domain, immunoglobulin domain (Ig), transmembrane domain (TM) and short cytoplasmic domain, (semaphorin) 4C 20353 ENSMUSG0….
Il1rl2|Il1rl2 153.74171 23.3046391 3.5649947 6.537075 0.0000000 0.0000000 Il1rl2 interleukin 1 receptor-like 2 107527 ENSMUSG0….
Fam126b|Fam126b 8348.97120 -0.8264253 0.2373035 -3.482567 0.0004966 0.0129644 Fam126b NA NA
Pard3bos2|Pard3bos2 48.12983 -23.7869821 3.6316680 -6.549878 0.0000000 0.0000000 Pard3bos2 par-3 family cell polarity regulator beta, opposite strand 2 105243940 ENSMUSG0….
Zdbf2|Zdbf2 3291.47772 -0.9550036 0.3014901 -3.167612 0.0015370 0.0316473 Zdbf2 zinc finger, DBF-type containing 2 73884 ENSMUSG0….
Klf7|Klf7 3114.99312 -0.6650797 0.2115799 -3.143398 0.0016700 0.0337915 Klf7 Kruppel-like factor 7 (ubiquitous) 93691 ENSMUSG0….
Map2|Map2 15113.61849 -0.6774677 0.1781094 -3.803660 0.0001426 0.0048256 Map2 microtubule-associated protein 2 17756 ENSMUSG0….
Arpc2|Arpc2 3532.07869 0.7956336 0.2237915 3.555245 0.0003776 0.0104712 Arpc2 actin related protein 2/3 complex, subunit 2 76709 ENSMUSG0….
Cnppd1|Cnppd1 1597.07545 0.9103458 0.3024970 3.009437 0.0026173 0.0475956 Cnppd1 cyclin Pas1/PHO80 domain containing 1 69171 ENSMUSG0….
Dnajb2|Dnajb2 1487.62355 1.1910669 0.3095799 3.847364 0.0001194 0.0041384 Dnajb2 DnaJ heat shock protein family (Hsp40) member B2 56812 ENSMUSG0….
Atg9a|Atg9a 8813.08524 2.0179228 0.4649912 4.339701 0.0000143 0.0006435 Atg9a autophagy related 9A 245860 ENSMUSG0….
Resp18|Resp18 56768.10866 -0.8522299 0.2302091 -3.701983 0.0002139 0.0066511 Resp18 regulated endocrine-specific protein 18 19711 ENSMUSG0….
Chpf|Chpf 11304.97381 2.0829954 0.3845636 5.416518 0.0000001 0.0000044 Chpf chondroitin polymerizing factor 74241 ENSMUSG0….
Tmem198|Tmem198 3311.19558 1.5407761 0.3702559 4.161382 0.0000316 0.0013071 Tmem198 transmembrane protein 198 319998 ENSMUSG0….
Slc4a3|Slc4a3 11951.48609 0.8940682 0.1933665 4.623698 0.0000038 0.0001938 Slc4a3 solute carrier family 4 (anion exchanger), member 3 20536 ENSMUSG0….

1.3 Visualize Shrinkage Estimations

Following subsetting, I want to see how the data looks based on effect size so I have run a shrinkage model using the “normal” method which normalizes counts based on effect size without major changes to data structure. I have included non-shrunk models as well for comparison.

1.3.1 Non-Shrunken L2FC

1.3.1.1 Craving vs FoodTrained

Version Author Date
dc04b15 avm27 2024-07-17

1.3.1.2 Maintenance vs FoodTrained

Version Author Date
dc04b15 avm27 2024-07-17

1.3.1.3 Craving vs Maintenance

Version Author Date
dc04b15 avm27 2024-07-17

1.3.2 Shrunken L2FC

1.3.2.1 Craving vs FoodTrained

Version Author Date
dc04b15 avm27 2024-07-17

1.3.2.2 Maintenance vs FoodTrained

Version Author Date
dc04b15 avm27 2024-07-17

1.3.2.3 Craving vs Maintenance

Version Author Date
dc04b15 avm27 2024-07-17

2 Looking at Differentially Expressed Genes

We can generate a matrix to identify where our differentially expressed genes lie between each group based on an alpha value of 0.05.

## Matrix Visualizattion
x <- vsDEGMatrix(
  data = dds, padj = 0.05, d.factor = "condition",
  type = "deseq", title = TRUE, legend = TRUE, grid = TRUE
)

Version Author Date
dc04b15 avm27 2024-07-17

3 Data Visualization

3.1 Normalized Counts and Plots of Genes of Interest

level_order <- c("FoodTrained", "Maintenance", "Craving") # this vector might be useful for other plots/analyses

normalized_counts <- counts(dds, normalized = TRUE)

3.1.1 Grin1

#### Grin1

d2 <- plotCounts(dds,
  gene = "Grin1|Grin1", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Grin1 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.1.2 Gabbr1

#### Gabbr1


d2 <- plotCounts(dds,
  gene = "Gabbr1|Gabbr1", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Gabbr1 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.1.3 Th

#### Th


d2 <- plotCounts(dds,
  gene = "Th|Th", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Th Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.1.4 Ehmt2

#### Ehmt2

d2 <- plotCounts(dds,
  gene = "Ehmt2|Ehmt2", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Ehmt2 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.1.5 Setdb1

#### Setdb1

d2 <- plotCounts(dds,
  gene = "Setdb1|Setdb1", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Setdb1 Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.1.6 Atf7ip

#### Atf7ip

d2 <- plotCounts(dds,
  gene = "Atf7ip|Atf7ip", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Atf7ip Expression",
    subtitle = "Normalized Gene Expression"
  )

a

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.1.7 Ash1l

#### Ash1l

d2 <- plotCounts(dds,
  gene = "Ash1l|Ash1l", intgroup = "condition",
  returnData = TRUE
)

a <- ggplot(d2, aes(x = factor(condition, level = level_order), y = count, color = condition)) +
  geom_point(position = position_jitter(w = 0.1, h = 0), size = 5) +
  scale_color_manual(values=c("#0000FF", "#FF8000", "#AD07E3")) +
  theme_bw() +
  xlab("Condition")+
  ylab("Normalized Counts")+
  ggtitle(
    label = "Ash1l Expression",
    subtitle = "Normalized Gene Expression"
  )

a

3.2 Principal Components Analysis (PCA)

For visualization I am converting all analyses to logarithmic. Following this I have completed a PCA for each analyses based on all genes together.

Version Author Date
bfc97c1 avm27 2024-07-17
1ee1d08 avm27 2024-07-17
41cc00d avm27 2024-07-17
741e5c2 avm27 2024-07-17
23fe19a avm27 2024-07-17
6284345 avm27 2024-07-17
5bc128a avm27 2024-07-17
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
cb49269 avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.3 Venn Diagrams

3.3.1 Up-Regulated Genes

myCol <- brewer.pal(3, "Pastel2")

results_sig_Table_up <- subset(results_sig, log2FoldChange > 0)
results_sig2_Table_up <- subset(results_sig2, log2FoldChange > 0)
results_sig3_Table_up <- subset(results_sig3, log2FoldChange > 0)

new_ressig1_up  <- rownames(results_sig_Table_up)
new_ressig2_up <- rownames(results_sig2_Table_up)
new_ressig3_up  <- rownames(results_sig3_Table_up)

# Generate Venn Diagram 
venn.plot.up <- venn.diagram(
  x = list(
    "Crv_vs_FT" = new_ressig1_up,
    "Maint_vs_FT" = new_ressig2_up,
    "Crv_vs_Maint" = new_ressig3_up),
    fill = myCol,
    cex = .8,
    fontface = "bold",
    fontfamily = "sans",
    cat.cex = 0.8,
    cat.fontface = "bold",
    cat.default.pos = "outer",
    cat.pos = c(0, 0, 0),
    cat.fontfamily = "sans",
    rotation = 1,
  category.names = c("Craving_vs_FoodTrained", "Maintenance_vs_FoodTrained", "Craving_vs_Maintenance"),
  filename = NULL,
  output = TRUE
)

# Save Venn Diagram
#pdf("plots/CIVSA_venn_diagram_Up.pdf")
#grid.draw(venn.plot.up)
#dev.off()

# Display Venn Diagram in RStudio
grid.draw(venn.plot.up)

Version Author Date
dc04b15 avm27 2024-07-17

3.3.2 Down-Regulated Genes

results_sig_Table_down <- subset(results_sig, log2FoldChange < 0)
results_sig2_Table_down <- subset(results_sig2, log2FoldChange < 0)
results_sig3_Table_down <- subset(results_sig3, log2FoldChange < 0)

new_ressig1_down  <- rownames(results_sig_Table_down)
new_ressig2_down  <- rownames(results_sig2_Table_down)
new_ressig3_down  <- rownames(results_sig3_Table_down)


# Generate Venn Diagram UP
venn.plot.down <- venn.diagram(
  x = list(
    "Crv_vs_FT" = new_ressig1_down,
    "Maint_vs_FT" = new_ressig2_down,
    "Crv_vs_Maint" = new_ressig3_down),
    fill = myCol,
    cex = .8,
    fontface = "bold",
    fontfamily = "sans",
    cat.cex = 0.8,
    cat.fontface = "bold",
    cat.default.pos = "outer",
    cat.pos = c(0, 0, 0),
    cat.fontfamily = "sans",
    rotation = 1,
  category.names = c("Craving_vs_FoodTrained", "Maintenance_vs_FoodTrained", "Craving_vs_Maintenance"),
  filename = NULL,
  output = TRUE
)

# Save Venn Diagram
#pdf("plots/CIVSA_venn_diagram_Down.pdf")
#grid.draw(venn.plot.down)
#dev.off()

# Display Venn Diagram in RStudio
grid.draw(venn.plot.down)

Version Author Date
dc04b15 avm27 2024-07-17

3.4 Heatmaps of Differentially Expressed Genes

ressig1 <- as.data.frame(results_sig)
ressig2 <- as.data.frame(results_sig2)
ressig3 <- as.data.frame(results_sig3)

ressig1$Comparison <- "Craving vs FoodTrained"
ressig2$Comparison <- "Maintenance vs FoodTrained"
ressig3$Comparison <- "Craving vs Maintenance"

ressig1_highlfc_pos <- ressig1[ressig1$log2FoldChange >= 1.5,]
ressig1_highlfc_neg <- ressig1[ressig1$log2FoldChange <= -1.5,]
ressig2_highlfc_pos <- ressig2[ressig2$log2FoldChange >= 1.5,]
ressig2_highlfc_neg <- ressig2[ressig2$log2FoldChange <= -1.5,]
ressig3_highlfc_pos <- ressig3[ressig3$log2FoldChange >= 1.5,]
ressig3_highlfc_neg <- ressig3[ressig3$log2FoldChange <= -1.5,]

temp_sigdegshm <- rbind(ressig1_highlfc_pos, ressig1_highlfc_neg)
temp_sigdegshm2 <- rbind(temp_sigdegshm, ressig2_highlfc_pos)
temp_sigdegshm3 <- rbind(temp_sigdegshm2, ressig2_highlfc_neg)
temp_sigdegshm4 <- rbind(temp_sigdegshm3, ressig3_highlfc_pos)
temp_sigdegshm5 <- rbind(temp_sigdegshm4, ressig3_highlfc_neg)

temp_sigdegshm6 <- temp_sigdegshm5[temp_sigdegshm5$lfcSE <=1,]

dup_rows <- duplicated(temp_sigdegshm6$symbol)
data_dedup <- temp_sigdegshm6[!dup_rows, ]

rows <- rownames(data_dedup)

select <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)#[1:1500]

mat <- assay(ddsMat_rlog)[select,]

data_mod <- mat[rownames(mat) %in% rows, ] 

matrownames_data_mod <- data.frame(rownames(data_mod))

matrownames_data_mod <- matrownames_data_mod %>% separate(rownames.data_mod., c("id", "symbol"),sep ="([|])")

samples <- as.tibble(read.csv("data/sample_info_CIVSA.csv", sep = ","))

annot_col <- samples %>%
  column_to_rownames('samplename') %>%
  as.data.frame()

ann_colors = list(
  condition = c(FoodTrained = "#0000FF", Maintenance = "#AD07E3", Craving = "#FF8000"))

3.4.1 Unsupervised Clustering

pheatmap::pheatmap(
  mat = data_mod,
  color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(155),
  scale = "row",
  annotation_colors = ann_colors,
  annotation_col = annot_col,
  annotation_legend = T,
  fontsize = 10, 
  show_colnames = T,
  cluster_rows = T,
  cluster_cols = T,
  show_rownames = F,
  clustering_distance_rows = "correlation",
  clustering_method = "ward.D",
  main = "Heatmap of Significant DEGs -- Unsupervised Clustering",
  cluster_row_slices = TRUE,
  show_row_dend = FALSE,
  cutree_rows = 4,
  cutree_cols = 3,
  treeheight_row = 10,
  treeheight_col = 10,
  #filename = "plots/CIVSA_Heatmap_DEGs_Final_Condensed_Unsupervised.pdf",
  legend_labels = c("FoodTrained", "Maintenance", "Craving"),
)

Version Author Date
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.4.2 Supervised Clustering

pheatmap::pheatmap(
  mat = data_mod,
  color = colorRampPalette(rev(brewer.pal(11, "RdYlBu")))(155),
  scale = "row",
  annotation_col = annot_col,
  annotation_colors = ann_colors,
  fontsize = 10,
  show_colnames = T,
  cluster_rows = T,
  cluster_cols = F,
  annotation_legend = T,
  legend = T,
  show_rownames = F,
  clustering_distance_rows = "correlation",
  clustering_method = "ward.D",
  main = "Heatmap of Significant DEGs -- Supervised Clustering",
  cluster_row_slices = TRUE,
  show_row_dend = FALSE,
  cutree_rows = 4,
  cutree_cols = 3,
  gaps_col = c(3, 6),
  treeheight_row = 10,
  treeheight_col = 10,
  #filename = "plots/CIVSA_Heatmap_DEGs_Final_Condensed_Supervised.pdf",
  legend_labels = c("FoodTrained", "Maintenance", "Craving"),
)

Version Author Date
12fbe39 avm27 2024-07-17
25a76da avm27 2024-07-17
21b70e3 avm27 2024-07-17
a1c090c avm27 2024-07-17
dc04b15 avm27 2024-07-17

3.5 Volcano Plots

3.5.1 Visualize Results: Volcano Plots (Non-Interactive)

Volcano plots are generated by gathering the FDR corrected p-values from each analysis. The adjusted p-values undergo a -log10 transformation to generate these FDR corrected values. Any rows containing “NA” are ommitted from analysis. Data points are colored based on increasing or decreasing value and plotted using ggplot2. Adjusted p-value cutoff is 1.3 following log transformation.(-log10 0.05 = 1.3)

Version Author Date
a1c090c avm27 2024-07-17

3.5.2 Visualize Results: Volcano Plots (Interactive)

Interactive Volcano plots are generated by Using the same method as above, but are plotted so that genes may be identified and log2FC can be seen from each gene. Plots can be used to customize which genes appear based on what you would like to see.

4 Pathway Analysis

4.1 Pathway Analysis GO

First, the background for the GO universe must be set. In order to do so we will subset the first results file “res” for any genes that map to entrez IDs and then we only take the unique ones. This leaves us with ~31,000 genes for the remaining background which we list as our midbrain dopamine neuronal expressed genes.

Pathway analysis is conducted using enrichGo and enrichKEGG based on a p-value of 0.05 and a q-value of 0.10 and including all significantly differentially expressed genes.

# Remove any genes that have greater than 1.5 LFC Standard Error
results_sig_entrez <- subset(results_sig, lfcSE <= 1.5)
results_sig_entrez2 <- subset(results_sig2, lfcSE <= 1.5)
results_sig_entrez3 <- subset(results_sig3, lfcSE <= 1.5)

# Remove any genes that have greater than adj.pvalue 0.05
results_sig_entrez <- subset(results_sig_entrez, padj <= 0.05)
results_sig_entrez2 <- subset(results_sig_entrez2, padj <= 0.05)
results_sig_entrez3 <- subset(results_sig_entrez3, padj <= 0.05)

# Remove any genes that do not have any entrez identifiers
results_sig_entrez <- subset(results_sig_entrez, is.na(entrez) == FALSE)
results_sig_entrez2 <- subset(results_sig_entrez2, is.na(entrez) == FALSE)
results_sig_entrez3 <- subset(results_sig_entrez3, is.na(entrez) == FALSE)

# Create a matrix of gene log2 fold changes
gene_matrix <- results_sig_entrez$log2FoldChange
gene_matrix2 <- results_sig_entrez2$log2FoldChange
gene_matrix3 <- results_sig_entrez3$log2FoldChange

# Add the entrezID's as names for each logFC entry
names(gene_matrix) <- results_sig_entrez$entrez
names(gene_matrix2) <- results_sig_entrez2$entrez
names(gene_matrix3) <- results_sig_entrez3$entrez


go <- enrichGO(
  gene = names(gene_matrix),
  OrgDb = "org.Mm.eg.db",
  readable = T,
  universe = gene_matrix_background,
  ont = "ALL",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.10
)

go2 <- enrichGO(
  gene = names(gene_matrix2),
  OrgDb = "org.Mm.eg.db",
  readable = T,
  universe = gene_matrix_background,
  ont = "ALL",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.10
)

go3 <- enrichGO(
  gene = names(gene_matrix3),
  OrgDb = "org.Mm.eg.db",
  readable = T,
  universe = gene_matrix_background,
  ont = "ALL",
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.10
)

4.1.1 Craving vs FoodTrained

kable(head(go, 10), caption = "Top 10 GO Pathways for Craving vs FoodTrained", format = "html", align = "l", booktabs = TRUE) %>%
  kable_styling(html_font = "Arial", font_size = 10, bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  scroll_box(width = "1000px", height = "500px")
Top 10 GO Pathways for Craving vs FoodTrained
ONTOLOGY ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
GO:0005344 MF GO:0005344 oxygen carrier activity 2/50 11/20962 0.0003025 0.0208720 0.0168759 Hba-a1/Hbb-bt 2
GO:0046974 MF GO:0046974 histone methyltransferase activity (H3-K9 specific) 2/50 11/20962 0.0003025 0.0208720 0.0168759 Ash1l/Setdb1 2
GO:0016209 MF GO:0016209 antioxidant activity 3/50 75/20962 0.0007641 0.0329752 0.0266618 Hba-a1/Nxn/Hbb-bt 3
GO:0008170 MF GO:0008170 N-methyltransferase activity 3/50 81/20962 0.0009558 0.0329752 0.0266618 Pemt/Ash1l/Setdb1 3
GO:0019825 MF GO:0019825 oxygen binding 2/50 23/20962 0.0013662 0.0377084 0.0304889 Hba-a1/Hbb-bt 2
GO:0051087 MF GO:0051087 chaperone binding 3/50 106/20962 0.0020727 0.0476730 0.0385457 Sdf2l1/Tfrc/H2-K1 3

4.1.2 Maintenance vs FoodTrained

kable(head(go2, 10), caption = "Top 10 GO Pathways for Maintenance vs FoodTrained", format = "html", align = "l", booktabs = TRUE) %>%
  kable_styling(html_font = "Arial", font_size = 10, bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  scroll_box(width = "1000px", height = "500px")
Top 10 GO Pathways for Maintenance vs FoodTrained
ONTOLOGY ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
GO:0031346 BP GO:0031346 positive regulation of cell projection organization 36/660 429/21359 1.00e-07 0.0002952 0.0002658 Arpc2/Obsl1/Cd24a/Palm/Cobl/Dlg4/Itga3/Mapt/Sphk1/Hnrnpk/Plxnb2/Twf1/Opa1/Ptk7/Ppp2r5b/Caprin1/Src/Mfn1/Ptprd/Wrap73/Ppp1r35/Fscn1/Rac1/Met/Megf8/Ap2a1/Ptpn5/Rgma/Apbb1/Cx3cl1/Katnb1/Ndrg4/Trim67/Sema7a/Map2k1/Golga4 36
GO:0010976 BP GO:0010976 positive regulation of neuron projection development 22/660 222/21359 1.60e-06 0.0036604 0.0032953 Cd24a/Cobl/Itga3/Mapt/Sphk1/Hnrnpk/Plxnb2/Twf1/Opa1/Ptk7/Ppp2r5b/Caprin1/Mfn1/Met/Ap2a1/Ptpn5/Rgma/Apbb1/Cx3cl1/Katnb1/Ndrg4/Trim67 22
GO:0098742 BP GO:0098742 cell-cell adhesion via plasma-membrane adhesion molecules 21/660 213/21359 3.10e-06 0.0045449 0.0040916 Obsl1/Cadm3/Cd24a/Cdh24/Pcdh8/Plxnb2/Igsf11/Ptprs/Pcdha4/Mbp/Cdh22/Fat4/Ptprd/Clstn3/Nectin2/Cadm4/Lrrc4b/Ric8a/Cx3cl1/Cdh13/Pcdh19 21
GO:0046777 BP GO:0046777 protein autophosphorylation 22/660 241/21359 6.30e-06 0.0056489 0.0050855 Clk1/Mapkapk2/Cdk12/Aatk/Calm1/Map3k1/Pim3/Acvr1b/Dyrk1a/Ddr1/Ppp2r5b/Map3k11/Stk39/Csnk2a1/Src/Jun/Epha8/Pink1/Atp13a2/Met/Calm3/Ulk3 22
GO:0071805 BP GO:0071805 potassium ion transmembrane transport 20/660 206/21359 6.60e-06 0.0056489 0.0050855 Calm1/Kcnk9/Kcnh8/Kcnk12/Slc12a2/Stk39/Slc12a5/Kcnq2/Hcn3/Kcna2/Slc24a2/Kcnab2/Kcnh2/Abcb8/Kcnip4/Pkd2/Calm3/Slc9a5/Nedd4/Kcnd1 20
GO:0050808 BP GO:0050808 synapse organization 34/660 490/21359 1.08e-05 0.0056489 0.0050855 Obsl1/Ngef/Palm/Dlg4/Cdk5r1/Itga3/Mapt/Rock2/Hnrnpk/Pcdh8/Ywhaz/Adgrb1/Plxnb2/Etv5/Opa1/Ptprs/Pcdhgc5/Sptbn2/Nrxn2/Caprin1/Mfn1/Bcan/Vcp/Ptprd/Ssh1/Setd5/Iqsec3/Clstn3/Lrrc4b/Ppfia3/Apbb1/Nedd4/Tpbg/Nlgn3 34
GO:0007156 BP GO:0007156 homophilic cell adhesion via plasma membrane adhesion molecules 14/660 114/21359 1.14e-05 0.0056489 0.0050855 Obsl1/Cadm3/Cdh24/Pcdh8/Plxnb2/Igsf11/Pcdha4/Cdh22/Fat4/Clstn3/Nectin2/Cadm4/Cdh13/Pcdh19 14
GO:0018107 BP GO:0018107 peptidyl-threonine phosphorylation 14/660 114/21359 1.14e-05 0.0056489 0.0050855 Clk1/Cdc42bpa/Dyrk2/Cdk5r1/Sphk1/Rock2/Dmtn/Acvr1b/Dyrk1a/Stk39/Csnk2a1/Dgkq/Met/Cadm4 14
GO:0032271 BP GO:0032271 regulation of protein polymerization 19/660 196/21359 1.14e-05 0.0056489 0.0050855 Arpc2/Camsap2/Inpp5j/Mapt/Hsp90aa1/Map3k1/Dmtn/Twf1/Dyrk1a/Snx9/Prkce/Ssh1/Arpc3/Rac1/Met/Mtpn/Eml2/Rps3/Cotl1 19
GO:0043254 BP GO:0043254 regulation of protein-containing complex assembly 30/660 410/21359 1.30e-05 0.0057749 0.0051989 Arpc2/Camsap2/Cd24a/Ncln/Inpp5j/Mapt/Hsp90aa1/Map3k1/Dmtn/Twf1/Wnt10b/Dyrk1a/Snx9/Prkce/Ppp2r5b/Src/Arhgef2/Vcp/Pink1/Ssh1/Arpc3/Stx1a/Fscn1/Rac1/Met/Mtpn/Eml2/Rps3/Prrt2/Cotl1 30

4.1.3 Craving vs Maintenance

kable(head(go3, 10), caption = "Top 10 GO Pathways for Craving vs Maintenance", format = "html", align = "l", booktabs = TRUE) %>%
  kable_styling(html_font = "Arial", font_size = 10, bootstrap_options = c("striped", "hover", "condensed", "responsive")) %>%
  scroll_box(width = "1000px", height = "500px")
Top 10 GO Pathways for Craving vs Maintenance
ONTOLOGY ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
GO:0046777 BP GO:0046777 protein autophosphorylation 22/604 241/21359 1.50e-06 0.0038779 0.0034854 Epha4/Mapkapk2/Mink1/Cdk12/Aatk/Calm1/Pim3/Acvr1b/Map3k12/Ddr1/Calm2/Ppp2r5b/Stk39/Src/Jun/Epha8/Pink1/Atp13a2/Met/Calm3/Hspa8/Ulk3 22
GO:0071805 BP GO:0071805 potassium ion transmembrane transport 20/604 206/21359 1.80e-06 0.0038779 0.0034854 Ywhae/Calm1/Kcnk9/Dlg1/Kcnk12/Calm2/Bin1/Stk39/Slc12a5/Hcn3/Kcna2/Kcnab2/Akap9/Kcnh2/Abcb8/Kcnip4/Calm3/Slc9a5/Nedd4/Kcnd1 20
GO:0035304 BP GO:0035304 regulation of protein dephosphorylation 13/604 96/21359 3.10e-06 0.0045037 0.0040478 Ywhae/Ppp1r9b/Ppp1r1b/Calm1/Ppp1r16a/Hsp90ab1/Calm2/Nsmf/Ptpa/Cry2/Pink1/Styxl1/Calm3 13
GO:0007041 BP GO:0007041 lysosomal transport 14/604 117/21359 5.70e-06 0.0062869 0.0056506 Plekhm1/Chmp7/Vps52/Bin1/Epg5/Vps51/Vps18/Snx16/Pink1/Atp13a2/Gak/Cacng7/Hspa8/Nedd4 14
GO:0035303 BP GO:0035303 regulation of dephosphorylation 15/604 136/21359 7.30e-06 0.0064833 0.0058271 Ywhae/Ppp1r9b/Ppp1r1b/Calm1/Ppp1r16a/Hsp90ab1/Calm2/Nsmf/Ptpa/Cry2/Src/Pink1/Styxl1/Ppp1r35/Calm3 15
GO:0031346 BP GO:0031346 positive regulation of cell projection organization 29/604 429/21359 1.53e-05 0.0112879 0.0101454 Arpc2/Epha4/Palm/Plekhm1/Rala/Twf1/Rrn3/Ppp2r5b/Grin1/Src/Stmn2/Wrap73/Fzd1/Styxl1/Ppp1r35/Fscn1/Met/Megf8/Ap2a1/Ptpn5/Rgma/Apbb1/Mob2/Crtc1/Cx3cl1/Katnb1/Ndrg4/Trim67/Nckipsd 29
GO:0016570 BP GO:0016570 histone modification 30/604 456/21359 1.84e-05 0.0116091 0.0104341 Elk4/Usp21/Ube2b/Srebf1/Asxl2/Hdac9/Bap1/Nipbl/Crebbp/Glyr1/Phf1/Ring1/Ubr2/Nfya/Dpy30/Cxxc1/Naa40/Ddb1/Phf20/Ash1l/Pygo2/Setdb1/Pink1/Ep400/Prmt8/Apbb1/Sin3b/Kmt2a/Setd2/Hdac8 30
GO:0006813 BP GO:0006813 potassium ion transport 20/604 243/21359 2.10e-05 0.0116091 0.0104341 Ywhae/Calm1/Kcnk9/Dlg1/Kcnk12/Calm2/Bin1/Stk39/Slc12a5/Hcn3/Kcna2/Kcnab2/Akap9/Kcnh2/Abcb8/Kcnip4/Calm3/Slc9a5/Nedd4/Kcnd1 20
GO:0048168 BP GO:0048168 regulation of neuronal synaptic plasticity 10/604 72/21359 3.36e-05 0.0148962 0.0133885 Nog/Grin2c/Syngr1/Nsmf/Grin1/Ncdn/Vgf/Rab8a/Jph3/Kmt2a 10
GO:0043666 BP GO:0043666 regulation of phosphoprotein phosphatase activity 9/604 58/21359 3.37e-05 0.0148962 0.0133885 Ppp1r9b/Ppp1r1b/Calm1/Hsp90ab1/Calm2/Ptpa/Cry2/Styxl1/Calm3 9

5 GO Pathways Visualization

5.1 GO DotPlot

Version Author Date
5bc128a avm27 2024-07-17

sessionInfo()
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] plotly_4.10.2                            
 [2] gridExtra_2.3                            
 [3] ggrepel_0.9.3                            
 [4] BiocParallel_1.32.6                      
 [5] conflicted_1.2.0                         
 [6] VennDiagram_1.7.3                        
 [7] futile.logger_1.4.3                      
 [8] regionReport_1.32.0                      
 [9] kableExtra_1.3.4                         
[10] vidger_1.18.0                            
[11] tables_0.9.17                            
[12] knitr_1.43                               
[13] Glimma_2.8.0                             
[14] Mus.musculus_1.3.1                       
[15] TxDb.Mmusculus.UCSC.mm10.knownGene_3.10.0
[16] OrganismDbi_1.40.0                       
[17] GenomicFeatures_1.50.4                   
[18] gplots_3.1.3                             
[19] lubridate_1.9.2                          
[20] forcats_1.0.0                            
[21] stringr_1.5.1                            
[22] purrr_1.0.2                              
[23] readr_2.1.4                              
[24] tibble_3.2.1                             
[25] tidyverse_2.0.0                          
[26] statmod_1.5.0                            
[27] tweeDEseqCountData_1.36.0                
[28] edgeR_3.40.2                             
[29] limma_3.54.2                             
[30] GOSemSim_2.24.0                          
[31] apeglm_1.20.0                            
[32] ggpubr_0.6.0                             
[33] SPIA_2.50.0                              
[34] KEGGgraph_1.58.3                         
[35] ggnewscale_0.4.9                         
[36] enrichplot_1.18.4                        
[37] tidyr_1.3.1                              
[38] WGCNA_1.72-1                             
[39] fastcluster_1.2.3                        
[40] dynamicTreeCut_1.63-1                    
[41] pathview_1.38.0                          
[42] gage_2.48.0                              
[43] dplyr_1.1.4                              
[44] topGO_2.50.0                             
[45] SparseM_1.81                             
[46] graph_1.76.0                             
[47] GO.db_3.16.0                             
[48] RColorBrewer_1.1-3                       
[49] genefilter_1.80.3                        
[50] pheatmap_1.0.12                          
[51] ggsci_3.0.0                              
[52] tximport_1.26.1                          
[53] org.Mm.eg.db_3.16.0                      
[54] AnnotationDbi_1.60.2                     
[55] DOSE_3.24.2                              
[56] ReactomePA_1.42.0                        
[57] biomaRt_2.54.1                           
[58] clusterProfiler_4.6.2                    
[59] ggplot2_3.5.0.9000                       
[60] DESeq2_1.38.3                            
[61] SummarizedExperiment_1.28.0              
[62] Biobase_2.58.0                           
[63] MatrixGenerics_1.10.0                    
[64] matrixStats_1.0.0                        
[65] GenomicRanges_1.50.2                     
[66] GenomeInfoDb_1.34.9                      
[67] IRanges_2.32.0                           
[68] S4Vectors_0.36.2                         
[69] BiocGenerics_0.44.0                      
[70] gprofiler2_0.2.2                         
[71] workflowr_1.7.0                          

loaded via a namespace (and not attached):
  [1] Hmisc_5.1-0              svglite_2.1.1            ps_1.7.5                
  [4] Rsamtools_2.14.0         foreach_1.5.2            rprojroot_2.0.3         
  [7] crayon_1.5.2             MASS_7.3-60              nlme_3.1-162            
 [10] backports_1.4.1          impute_1.72.3            rlang_1.1.3             
 [13] XVector_0.38.0           HDO.db_0.99.1            callr_3.7.3             
 [16] filelock_1.0.2           rjson_0.2.21             bit64_4.0.5             
 [19] glue_1.7.0               rngtools_1.5.2           parallel_4.2.1          
 [22] processx_3.8.2           DEFormats_1.26.0         tidyselect_1.2.1        
 [25] XML_3.99-0.14            GenomicAlignments_1.34.1 xtable_1.8-4            
 [28] magrittr_2.0.3           evaluate_0.21            bibtex_0.5.1            
 [31] cli_3.6.2                zlibbioc_1.44.0          doRNG_1.8.6             
 [34] rstudioapi_0.14          whisker_0.4.1            bslib_0.5.0             
 [37] rpart_4.1.19             derfinderHelper_1.32.0   fastmatch_1.1-3         
 [40] BiocStyle_2.26.0         lambda.r_1.2.4           treeio_1.22.0           
 [43] xfun_0.39                gson_0.1.0               cluster_2.1.4           
 [46] caTools_1.18.2           tidygraph_1.2.3          KEGGREST_1.38.0         
 [49] ape_5.7-1                Biostrings_2.66.0        png_0.1-8               
 [52] reshape_0.8.9            withr_3.0.0              bitops_1.0-7            
 [55] ggforce_0.4.1            RBGL_1.74.0              plyr_1.8.8              
 [58] coda_0.19-4              bumphunter_1.40.0        pillar_1.9.0            
 [61] cachem_1.0.8             fs_1.6.2                 graphite_1.44.0         
 [64] ellipsis_0.3.2           vctrs_0.6.5              generics_0.1.3          
 [67] tools_4.2.1              foreign_0.8-84           munsell_0.5.1           
 [70] tweenr_2.0.2             fgsea_1.24.0             DelayedArray_0.24.0     
 [73] fastmap_1.1.1            compiler_4.2.1           abind_1.4-5             
 [76] httpuv_1.6.11            rtracklayer_1.58.0       GenomeInfoDbData_1.2.9  
 [79] lattice_0.21-8           utf8_1.2.4               later_1.3.1             
 [82] BiocFileCache_2.6.1      jsonlite_1.8.7           GGally_2.1.2            
 [85] scales_1.3.0             tidytree_0.4.2           carData_3.0-5           
 [88] lazyeval_0.2.2           promises_1.2.0.1         car_3.1-2               
 [91] doParallel_1.0.17        checkmate_2.2.0          rmarkdown_2.23          
 [94] cowplot_1.1.1            webshot_0.5.5            downloader_0.4          
 [97] BSgenome_1.66.3          igraph_1.5.0             survival_3.5-5          
[100] numDeriv_2016.8-1.1      yaml_2.3.7               systemfonts_1.0.4       
[103] htmltools_0.5.5          memoise_2.0.1            VariantAnnotation_1.44.1
[106] BiocIO_1.8.0             locfit_1.5-9.8           graphlayouts_1.0.0      
[109] viridisLite_0.4.2        digest_0.6.32            rappdirs_0.3.3          
[112] knitrBootstrap_1.0.2     futile.options_1.0.1     emdbook_1.3.13          
[115] RSQLite_2.3.1            yulab.utils_0.0.6        derfinder_1.32.0        
[118] data.table_1.14.8        blob_1.2.4               preprocessCore_1.60.2   
[121] labeling_0.4.3           splines_4.2.1            Formula_1.2-5           
[124] RCurl_1.98-1.12          broom_1.0.5              hms_1.1.3               
[127] colorspace_2.1-0         base64enc_0.1-3          BiocManager_1.30.21     
[130] aplot_0.1.10             nnet_7.3-19              sass_0.4.6              
[133] Rcpp_1.0.10              mvtnorm_1.2-2            fansi_1.0.6             
[136] tzdb_0.4.0               R6_2.5.1                 lifecycle_1.0.4         
[139] formatR_1.14             curl_5.0.1               ggsignif_0.6.4          
[142] jquerylib_0.1.4          Matrix_1.5-4.1           qvalue_2.30.0           
[145] org.Hs.eg.db_3.16.0      iterators_1.0.14         RefManageR_1.4.0        
[148] htmlwidgets_1.6.2        markdown_1.7             polyclip_1.10-4         
[151] crosstalk_1.2.0          shadowtext_0.1.2         timechange_0.2.0        
[154] gridGraphics_0.5-1       reactome.db_1.82.0       rvest_1.0.3             
[157] htmlTable_2.4.1          patchwork_1.2.0.9000     bdsmatrix_1.3-6         
[160] codetools_0.2-19         gtools_3.9.4             getPass_0.2-2           
[163] prettyunits_1.1.1        dbplyr_2.3.2             gtable_0.3.5            
[166] DBI_1.1.3                git2r_0.32.0             ggfun_0.1.1             
[169] httr_1.4.6               highr_0.10               KernSmooth_2.23-21      
[172] stringi_1.8.4            progress_1.2.2           reshape2_1.4.4          
[175] farver_2.1.2             annotate_1.76.0          viridis_0.6.3           
[178] Rgraphviz_2.42.0         ggtree_3.6.2             xml2_1.3.4              
[181] bbmle_1.0.25             restfulr_0.0.15          geneplotter_1.76.0      
[184] ggplotify_0.1.1          bit_4.0.5                scatterpie_0.2.1        
[187] ggraph_2.1.0             pkgconfig_2.0.3          rstatix_0.7.2           
[190] GenomicFiles_1.34.0